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Abstract 

We introduce a numerical method for extracting minimal geodesics along the group of 
volume preserving maps, equipped with the metric, which as observed by Arnold |Arn66| 
solve Euler’s equations of inviscid incompressible fluids. The method relies on the general¬ 
ized polar decomposition of Brenier |Bre m, numerically implemented through semi-discrete 
optimal transport. It is robust enough to extract non-classical, multi-valued solutions of Eu¬ 
ler’s equations, for which the flow dimension is higher than the domain dimension, a striking 
and unavoidable consequence of this model |Shn94) . Our convergence results encompass this 
generalized model, and our numerical experiments illustrate it for the first time in two space 
dimensions. 


1 Introduction 

The motion of an inviscid incompressible fluid, moving in a compact domain X C is described 
by Euler’s |Eul65| equations 

dtv + {v • V)u = —Vp divu = 0, (1) 

coupled with the impervious boundary condition u • n = 0 on SO. Here v denotes the fluid 
velocity, and p the pressure acts as a Lagrange multiplier for the incompressibility constraint. 
In Lagrangian coordinates, Euler equations 0 yield the geodesic equations along the group 
SDiff volume preserving diffeomorphisms of X, equipped with the metric |Arn66| . Consider 
an inviscid incompressible fluid flowing during the time interval [0,1], and a map 5* : X ^ X 
giving the final position s*(x) of each fluid particle initially at position x G X. In this paper, we 
discretize and numerically investigate a natural approach to reconstruct the intermediate fluid 
states: look for a minimizing geodesic joining the initial configuration — Id to the final one 5* 

minimize f \\s{t)\\‘^dt^ subject to s(0) = 5*, s(l) = s*, and Vt G [0,1], s{t) G S. (2) 

Jo 

We denoted by S C L^(X, M^) the space of maps preserving the Lebesgue measure on X, which 
in dimension d > 2 is the closure of SDiff. Despite this first relaxation, note that the optimized 
functional in 0 does not penalize the spatial derivatives of s, whereas the constraint involves 
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Figure 1: The motion of inviscid incompressible fluids admits three formulations, either (I) 
Eulerian based on the local speed u : [0,1] x X ^ (II) Lagrangian based on diffeomorphisms 
5(t, •) which integrate the speed: dts{t^x) — or (III) relaxed as a superposition of 

individual particles paths cj G ri, weighted by a measure /i. 


the Jacobian of s. The study of © thus requires non-standard variational techniques, reviewed 
in |FD12| . 

In dimension d > 3, the optimization problem © needs not have a minimizer in 5 G 
|Shn94| . and minimizing sequences {sn)neN instead display oscillations reminis¬ 
cent of an homogeneization phenomenon. A second relaxation is required, based on generalized 
flows [Bre89| which allow particles to split and their paths to cross. This surprising behavior is 
an unavoidable counterpart of the lack of viscosity in Euler’s equations, which amounts to an 
infinite Reynolds number. Generalized flows are also relevant in dimension d G {1,2} if the un¬ 
derlying physical model actually involves a three dimensional domain X x in which one 

neglects the fluid acceleration in the extra dimensions |Bre08| . Consider the space of continuous 
paths (of fluid particles) 

:=C°([0,1],X). 

Let et{uj) := uj{t) be the evaluation map at time t G [0,1], so that (eo, ei){(jj) = (cc;(0),Cc;(l)). Let 
also Leb denote the Lebesgue measure restricted to the domain X, normalized for unit mass, 
and let /#/i denote the push-forward of a measure /i by a measurable map /. The geodesic 
distance Q admits a convex relaxation, linearizing both the objective and the constraints, and 
for which the existence of a minimizer is guaranteed. It is posed on probability measures on fl, 
called generalized flows 


d^(5x,,s*):= min / subject to < 

/iGProb(O) Jq 


— Id 

(eo, ei)#/i = (s*, s*)# Leb, 
yt G [0,1], et = Leb. 


(3) 


Note that the path action A : ^ ^ M+ U {-hoo}, although unbounded, is lower semi-continuous. 
The first constraint (eo,ei)#/i = (5*,s*)#Leb expresses that moving fluid particles from s^{x) 
to 5*(x) for all X G X, or from the origin (x;(0) to the end of the paths cj G 12 as weighted 
by /i, yields equivalent transport plans. The second constraint e^ #/i = Leb states that the path 
positions as weighted by /i, equidistribute on X at each time t G [0,1], which amounts to 
incompressibility. A classical flow s G iL^([0,1],S) can be regarded as a generalized flow, with 
paths t i-G 5(t,x), weighted by the Lebesgue measure on x G X. Our discretization truly solves 
©, rather than Q, and convergence is established in this relaxed setting. 

The incompressibility constraint in ©, © and ©, gives rise to a Lagrange multiplier, 
the pressure, which is the unique maximizer to a concave optimization problem dual to Q, 
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Figure 2: The geodesic distance along the “manifold” S of volume preserving maps, 

represented as a blue curve, is estimated as the length of a chain {mpJ^Q in the linear subspace 
Mat, represented as a black line, plus penalizations for the boundary values and the distance from 
the chain elements to S. 


see |Bre 93]. The primal Q may in contrast have several solutions, up to the notable ex¬ 
ception |BFS09| of smooth flows in dimension d = 1. The pressure is a classical function 
p ^ Bv(x)), see [AF07) (which requires the technical assumption that X is a d- 

dimensional torus). This regularity is sufficient to show that any solution s to 0 (resp. /i-almost 
any path cj, for any solution /i to Q) satisfies 

dus{t,x) = -Vp(t, s(t,x)), resp. u{t) = (4) 

In other words, fluid particles move by inertia, only deflected by the force of pressure. Assume 
that the pressure hessian is sufficiently small, precisely that 

Vt G [0,1], Vx G X, x) -< tP Id (5) 

in the sense of symmetric matrices. Then using the path dynamics equation Q Brenier |Bre89| 
showed that the relaxed problem Q admits a unique minimizer p G Prob(i2), which is deter¬ 
ministic: in other words associated to a, possibly non-smooth but otherwise classical, minimizer 
seHH[0,Tl^) of©. Inequality ^ is sharp, and several families of examples are known for 
which uniqueness and/or determinism are lost precisely when the threshold © is passed. We 
present ^the first numerical illustration of this phenomenon. 


1.1 Numerical scheme and main results 

We introduce a new discretization for the relaxation of the shortest path formulation Q of 
Euler equations 0. Our approach is numerically tractable in dimension 2, and is the first to 
illustrate the transition between classical and generalized solutions occurring at the threshold 
0 on the pressure regularity. 

For that purpose we need to introduce some notation. Let M := L^(X, M^), and let S C M 
be the collection of maps preserving the restriction to X of the Lebesgue measure, denoted by 
Leb and normalized to have mass 1. For each X G N let Vn be a partition of X into N regions 
of equal area 1/X, diameter < CxN~d^ and let Mw C M be the X-dimensional subspace 
of functions which are piecewise constant on this partition. Given G S, discretization 

parameters T, X G N, and a penalization factor A » 1, we solve 


£{T,N,X) 


mSM^ 0 <i<r 


—5*^+ k«T —S P+ inf 

se§ 

i<i<r 



( 6 ) 

In all this paper, || • || stands for the L^(X, norm, and | • | for the euclidean norm on M^, for 
any n G N. Comparing this with Q, we recognize the standard discretization of the length of the 
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discrete path (mo, • • • , m^), as well as an implementation by penalization of the boundary value 
constraints and of the incompressibility constraints. The optimization of Q, seen as a function 
of m G is diU N{T -\- 1)d-dimensional smooth optimization problem. A quasi-Newton 

method gave convincing results, see il despite the non-convexity of the functional which forbids 
to guarantee that its global minimum is numerically found. 

Before entering the analysis of (§, let us emphasize that the inner-subproblems, the projec¬ 
tion of each rrii onto the set S of measure preserving maps, are numerically tractable thanks to 
two main ingredients: Brenier’s polar factorization [Bre91| . and semi-discrete optimal transport. 
The former states that the distance from any given m G M to the set S, is the cost of the 
transport plan needed to equidistribute on X the image measure of m 

inf \\m — s\\‘^ = W 2 {m# Leb, Leb), (7) 


where W 2 is the Wasserstein distance for the quadratic transport cost. If m G Mw, then m# Leb 
is the sum of N Dirac measures of mass 1/N, located at the N values of the piecewise constant 
map m on the partition Vn- Semi-discrete optimal transport |AHA98( IMerlfLevldl is a numer¬ 
ical method for computing Q, and more generally the Wasserstein distance between a discrete 
measure r] = and an absolutely continuous measure u = pLeb, with a (typically) 

piecewise linear density p. It is based on Kantorovitch duality 

0) = sup fv+ gi', where Vy € X, g{y) = inf \y - xf - f{x), (8) 

feL^(v)Jx Jx 

= sup V Vjfj+ 9^^ whereVye V ^(y)^ min (9) 

where Q is obtained from Q by setting fj = f{xj). Importantly, the conjugate y in Q is 
piecewise quadratic on a partition of X, called the Laguerre Diagram of the sites xj with weights 
/j, that is constructible through computational geometry software [cpij . The X-dimensional 
concave maximization problem (§, which is unconstrained and twice continuously differentiable, 
is efficiently solved via Newton or quasi-Newton methods. Semi-discrete optimal transport has 
become a reliable and efficient building block for PDE discretizations |BCM01^ . 

A second interpretation of the optimization problem (§, closer to Q , involves a generalized 
ffow pi G Prob(D) supported on N trajectories, each piecewise linear with direction changes at 
times {0,1/T, ••• ^TjT}. Let m — (m^)^Q G and for each 1 < j < X let mi be the 

constant value of on the j-th region of the partition Vn of X. For each 1 < j < X let 
ujj G D be the piecewise linear path with value m] at time i/T, for all 0 < i < T (see Figure]^. 
Finally let pi G Prob(D) be the discrete probability measure equidistributed on the set of paths 
{ujj] 1 < j < X}. Then Q rewrites in a form close to ^ 


f A{u)dpi{u) ^ \{ f \mQ{x) - s^{x)\‘^-\-\mT{x) — s''{x)\‘^dx-\- ^ lT|(et#/i, Leb) ]. (10) 
Jq \ Jx / 


l<i<T 

t=i/T 


Indeed, the first energy term satisfies 


f A{u)dfi{u) = ^ 


17 7 12 

mP 1 — mp = 


^<j<N 


z+l ll^^+l 

0<i<T 0<7<T 

^<3<N 


mi 
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Figure 3: (Left) A partition Vn cuts the domain X into N region of equal area and roughly 
isotropic shape. (Right) To a sequence (m^)^Q G one can associate N piecewise linear 

paths by interpolating the map values at the times {0,1/T, • • • /T} for each region 

of the partition Vn- 


The penalized integral term in ( [To| ) equals 


\Tno — 5 * 


+ II^T — 5*IP from Q. It is the cost of 


the transport plan on mapping {mQ{x)^mT{x)) to (sx,(x), s*(x)) for all x G X, which sends 
(eo,ei)#/i = (mo,mT)#Leb onto (s*,s*)#Leb, and thus enforces the proximity of these two 
couplings on X^ as required in Q. The other penalized terms VF|(et #/i, Leb) account for the 
incompressibility of /i at time t = i/T, l<i<T, and by Q are equal to inf^^g ||m^ — sjp from 

Summarizing, the geodesic formulation of Euler equations Q has a rather surprising re¬ 
laxation (§, looking a-priori unphysical: fluid particles may split and cross. Yet the natural 
discretizations and ( [To| ) of these two formulations are actually identical The classical and 
generalized interpretations are also at the heart of our main result. 

Theorem 1.1. Let G let T^N G N; and A > 0. The relaxed geodesic distance © and 

the discretized minimum Q satisfy 

£{T, N, A) < d{s^, s*) + OiTh%X), 

• (Classical estimate) with Ln = N~d^ if the classical geodesic distance © equals the relaxed 
distance ©, and admits a minimizer with regularity s G L^([0, l],iL^(X)). 

• (Relaxed estimate) with Ln = (resp. N~ 2 ^/\nN if d — \), if the pressure field 

gradient Vp is Lipschitz on [0,1] x X, and the boundary data 5“^, s* are Lipschitz on X. 

Recall that the classical Q and relaxed Q distances are automatically equal in dimension 
d > 3, and that the pressure held gradient Vp is uniquely determined by the boundary values 


The decay rate Ln = N ^ in Theorem 


1.1 


is actually tied to the dimension T^quant(M) of the 


generalized flow fi G Prob(ri) minimizing ([3), see Deflnition 3.1 The flow associated to a classical 
solution has dimension d, since the particle trajectories are determined by their initial position 
X G X C The trajectories of a generalized flow obey a second order ordinary differential 
equation 0 and are thus determined by their initial position and velocity (x, u) G X x C 
provided Cauchy-Lipschitz’s theorem applies. The generalized ffow dimension is thus 2d in the 
worst case, but intermediate dimensions d < D < 2d are also common, see 0 

Theorem LA does not tell how to choose the constraint penalization parameter A. The next 
proposition shows that the quantity S'{T^ X, A) := (1 -h4T/A)^^(r, X, A) arises naturally in error 
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( 11 ) 


estimates, which suggests to choose A = = Nd so that 

S\T, N, A) = s*) + 0{T/X + Th%X) = d\s^, 5 *) + 0{TN-^). 

Proposition 1.2. Let m = (m^)^Q G be a minimizer of (§. 

• (Classieal eonstruetion) Let (5^)^q be the ehain of ineompressible maps defined by: sq — 

St = and Si is a projeetion of rui onto S for all 1 < i <T. Then 

T \\si+i-sif <£'{T,N,X). 

0<i<T 

• (Relaxed eonstruetion) Let ji G Prob(fi) be the generalized flow built from (m^)^Q as in 
( [Tol ) . Then 

^ W^{et #//, Leb) dt < N, A). 

^5 a result, let (A^r, Ar)TGN be sueh that £\T,Nt^Xt) d{s^,s'') as T ^ 00. Then a 
subsequenee of the assoeiated flows (/iT)TGN weak-"^ eonverges to a minimizer of (§. 


Outline. Theorem O is established p.l| (Classical estimate) and ^ (Relaxed estimate). 
Proposition |1.2| is proved §2.2[ Numerical experiments are presented Q 


Remark 1.3 (Monge-Ampere gravitation). Some models of reeonstruetion of the early universe 
JBrell]/ involve aetions of a form elosely related to our diserete energy funetional (§, for the 
parameter value A = 2; 


^ (^ll™-(^)f+ inf ||m(t)-sf 


dt. 


2 Classical analysis 


We establish Theorem O (Classical estimate) in ^2.1| and prove Proposition |0| in ^2.2[ The 
optimization parameters (T, N, A, 5 *, 5 *) are fixed in this section. 


2.1 Upper estimate of the discretized energy 

Following the assumption of Theorem O (Classical estimate), we consider a minimizer of the 
shortest path problem and assume that it has regularity s G L^([0,1], i7^(X)). Define 
Si := s{i/T) for all 0 < i < T, and note that sq = s^, st = 5*. Let also rui := Fjsf^Si), for all 
0 < i < T, where Pw : M ^ Mw denotes the orthogonal projection. We denote by h^ := N~d 
the discretization scale, and recall that each region of the partition Vn of ^ has area 1/N and 
diameter < C^hN- 

Let s^ denote the mean of Si on the region P of the partition Vn, for all 0 < i < T. Then 


Si 




s[fdx<Csh{CrhN)‘^ ^ / \Vsi{x)fdx = Ch%\\'Vsi\\‘^, 
PeVN 

( 12 ) 


where the Sobolev inequality constant C^h only depends on the dimension, and C Cs\^C^. 
Recall that, in all this paper, || • || stands for the L^(X,M^) norm, and | • | for the euclidean norm 
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Figure 4: Theorem 

^ onto the finite dimensional space M^v, a procedure symmetric the projection of 

G onto S involved in the discrete energy optimization Q, see Figure 


1.1 


(classical estimate) is based projecting the measure preserving maps 


on for any integer n > 1. The map Pjv is 1-Lipschitz, as the orthogonal projection onto the 
convex set Mjv- Hence for any 0 < i < T 





s{t) W^dt. 


Summing (12) and ([T^ over 0 < i < T we obtain 


£{T,N,\) <T y] ||mi+i - + A \\mi-Si 

0<i<T 0<i<T 

— T/ /■ P(t)|pc?t + A y] Ch%\\Vsi 

0<i<T'' T 0<i<T 

<cf{s^,s*) + C'Th%X, 


(13) 


where C' — C'||5||2^ooQo,i],ifi(x))5 which concludes the proof. 

2.2 Length of a chain of incompressible maps 


Proposition |1.2| (Classical construction) immediately follows from Lemma 2.1 below, which is 
general and could be used to approximate geodesics on any manifold S embedded in a Hilbert 
space M, internally approximated by subspaces Mw- It relies on a the following identity, valid 
for any elements a, 6 of a Hilbert space, and any 6: > 0: 

{l + e)-^\\a + bf < \\af + e-^\\bf. (14) 

Indeed subtracting the LHS to the RHS of (14) we obtain (1 + e)~^\\e^a — e~ 2 b\\‘^ > 0. 
Lemma 2.1. For any T G N*; any penalization A > 0^ and any (m, s) G (M x S)^+^ one has 


T \\si+i-Sif <il + 4:T/X) 

0<i<T 


T y] ||mj+i-m*p + A \\mi - Si 
Q<i<T 0 <i<r 


(15) 


Proof. Let 0 < i < T. Choosing a := 5^+1 — 5^, and b := — a, we obtain 

(1 + £)~^\\mi+i - < pi+i - Sip + e“^|l(si - m*) - (sj+i - mi+i)p 

< Pi+i - Sip + 2e“p||si - mip + ||si+i - mi+ip). 
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Summing over 0 < i < T yields 


(l + e) ^ y] ||mi+i - mill^ < ||sj+i - + 2£ ^ ^ ai\\si-mi\\‘^, 

0<i<T 0<i<T 0<i<T 

with ao = aT — o^i — ‘2 otherwise. Choosing = 4T/A concludes the proof. □ 

The second point of Proposition |l.2| is based on ( [l4| ) as well. Indeed, let (m^)^Q be minimizers 
of Q, let 0 < i < T, and let t = {i + a)/T with 0 < a < 1. Then for any 6: > 0 

Leb) = inf ||(1 ——5 

Integrating over t G [0,1], using that either a < 1/2 or 1 — cr < 1/2, and choosing 6: = 4T we 
obtain as announced 

f W 2 (et#M,Leb)dt < ^ T||mi+i - inf ||mi - s||A = ^5'(T,iV, A). 

Finally, the convergence claim for the minimizing chain (/iT)rGN results from classical arguments, 
(i) The weak-* lower semi-continuity of the energy /i i-G A{(jj)diJi{(jj) on Prob(i2), which follows 
from the lower semi-continuity of the action ^ ^ U {oo}. (ii) The weak-* sequential 

compactness of {/i G Prob(i2); A(cj)dju(u;) < K} for any constant K, see |Bre93| . (hi) The 
weak-* continuity of /i i-G VFI((cq, ei)#/i, 5*)# Leb), a quantity bounded for fiT by ||mo — 

-h Wttit — s*|p < E{T^ Ntj Xt)/Xt ^ 0 as T ^ oc. (iv) The weak-* lower semi-continuity 
of /i i-G VFI(et #/i, Leb) dt, which follows from Fatou’s lemma and the continuity of /i i-G 
W^iet #^7 Leb) for any t G [0,1]. 


< (1+6:) ||o:(m^+i 


m^)|p +6: Mnf I 

seS 


rrii 


3 Relaxed analysis 

We prove Theorem 1 1.1 1 (Relaxed estimate), using a quantization of the generalized flow minimiz¬ 
ing the relaxed geodesic distance (§. This quantization is a counterpart of the partition Vn of 
the domain (X, Leb) used for the classical estimate ^2.1, which amounts to quantize the initial 
positions of the fluid particles. Let Sx denote the Dirac probability measure concentrated at a 
point X. 


Definition 3.1. Let IH 6e a metric space, let ji be a probability measure on and let F C M. 
For all N > 1 denote, with W 2 the Wasserstein distance for the quadratic transportation cost 

hNU)-.^ inf 1^2 (ft, f tnU) inf min|r>0; TC M B{ui,r)\. 


The quantization dimension of fi, and the box dimension of F^ are defined by 

In X In X 

L^quant(/^) := hmsup . , . L)box(r) := limsup —-—. 

w^oo -ln/iw(M) N^oo -lnrw(r) 

The decay rate of /iw is directly involved in the announced result Theorem |1.1[ We estimate 
it using an elementary result of quantization theory, and refer to |GG92| for more details on 
this rich subject. Note that the (upper) box dimension Dbox is a variant of the Haussdorff 
dimension, in which the set of interest if covered by balls of equal radius. Box and Haussdorff 











dimension coincide for compact manifolds, but differ in general. For instance, all countable sets 
have Haussdorff dimension zero, whereas one can check that 


^box(([0,i]nQ)'^) 


d. 



1 

2 


Proposition 3.2. Let M be a metrie spaee, and let ji G Prob(IHI) he supported on a set T. Then 
^quant(M) ^ max{2, D\^q^{T)}. More preeisely for any D > 0, one has as N ^ oo 



(16) 


Proof Let G N be fixed. For each 1 < i < let C IH be a set of i points such that 
r C 2r^), with := n(r). We construct a sequence of points oji G IH, and an 

increasing sequence of measures pi supported on F and of mass i/N^ inductively starting with 
i — N and finishing with i — 1. Initialization: p. 

Induction: for each 1 < i < A^, we construct oji and pi-i in terms of pi. Let indeed uji G Mi 
be such that Bi B{(jJi^2ri) satisfies Pi{Bi) > 1/N. Such a point exists since \pi\ = i/N^ 
#{Mi) = i, and supp(p^) C F. Then let pi-i := pi — so that pi — pi-i is a non-negative 

measure of mass ^ supported on Bi. One has 



l<i<N ^ l<i<N ^ ^ l<i<N 


l<i<N 


l<i<N 


The comparison (16) of the decay rates of h^ip) and rw(r) immediately follows. Finally the 


comparison of the dimensions follows from ( [T^ . 


□ 


We now specialize the choice of /i, F and IH. Let p G Prob(i2) be a generalized flow minimizing 
the relaxed geodesic distance This measure is concentrated on the set F of paths obeying 
Newton’s second law of motion 


T:={ue C\[0, 1],X); Vt G [0,1], u{t) = 


where the pressure gradient Vp : [0,1] x X ^ is assumed, following the assumptions of 
Theorem [HI to have Lipschitz regularity. We regard F as embedded in the Hilbert space 
IH := iL^([0,1], M^), which plays a natural role in the problem of interest (§ and is equipped 
with the norm 



Note that IH continuously embeds in (7^(12,M^), hence the evaluation maps e^ : IH ^ are 
continuous with a common Lipschitz constant denoted Cq. 

Lemma 3.3. The set F is eompaet. Furthermore the map F ^ X x : cj i-G ((x;(0),d;(0)) is 
bijeetive and bi-Lipsehitz onto its image. 

Proof. The result follows from Cauchy-Lipschitz’s theorem for ordinary differential equations, 
and the compactness of X. □ 
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The image of the generalized flow fi by the map of Lemma |3.3| namely initial position 
and speed, is often called a minimal measure |BFS09| . Since there is no ambiguity, we denote 
hN := /iw(m)- The constants c, C, C' appearing in the estimates below only depend on the 
dimension d. 


Corollary 3.4. One has hjsf — 0{N 2 d) (resp. 0{N ^y/luN) if d— 1.) 


Proof. By Lemma 3^, the set L is in bi-Lipschitz bijection with a compact set K C Hence 
^w(r) < CrN{K) < C'N~^, and the upper estimate follows from (16). □ 


The quantization scale h]\f is also bounded below, and is minimal for classical solutions. 

Lemma 3.5. There exists c > 0 sueh that hjsf > cN~d for all N > 0. If the generalized flow 
fi in faet represents a elassieal solution s to Euler^s equations, and Vs is bounded on [0,1] x X, 
then this lower estimate is sharp: h^ — 0{N~df 

Proof. Since X is a d-dimensional domain, there exists c > 0 such that lT2(Leb, z/w) ^ cN~d 
for any measure vjsi supported at N points of M^. (Recall that, in this paper, Leb denotes the 
Lebesgue measure restricted to the set X, and normalized for unit mass.) The first point follows: 
for any measure /iw supported at N points of IH 

cN~^ < lT2(Leb,eo#)niv) = bL2(eo < CeW2{ia,iaN) = Ce/iw- 


Second point: for each x G X, let t s{t,x). Then $ : (X, Leb) ^ (r,/i) : x i-G Ux is 
measure preserving and Lipschitz, with regularity constant denoted C^. Let z/w be a discrete 
probability measure, with one Dirac mass of weight 1/N in each region of the partition Vjy. 
Since these regions have diameter < C^N~d^ we conclude that 

hN < fT2(/i,$#z/iv) = lT2(^#Leb,$#z/w) < C<^lT2(Leb, z/^) < □ 


In the rest of this section, we fix the integer N and allow ourselves a slight abuse of notation: 
elements ujj, Pj^ Pj^ • • • indexed by 1 < j < X do implicitly depend on X, although that second 
index , P^, is omitted for readability. 


Lemma 3.6. The infimum defining h^ is attained, see Definition \3.1 
^ probability measures on T sueh that 


As a result there exists 


N Tj 

^<j<N 


h% — 




(jJ — CJn 


idpjH 


(17) 


Furthermore, ujj is the baryeenter of pj for eaeh 1 < j < X. 

Proof. Let be a candidate quantization, and let tt be the transport plan associated to 

^ujj^P)’ Then the measures pj : A \-^ N 7r{{xj} x A), 1 < j < X, are probabilities 
which average to /i, and the transport cost is the RHS of ( [T7| . The quantization energy, i.e. 
the squared Wasserstein distance, is decreased by replacing ojj with the baryeenter bj of pj, 
1 < J < X, by the amount ^ ^j=i \^j Hence cjj = bj for all 1 < j < X if the quantization 

is optimal. Note also that the baryeenter of pj belongs to G := Hull(r) by construction. 

Since T is a compact subset of a Hilbert space, the convex hull closure G is also compact, for 
the strong topology induced by || • ||e- The quantization energy i-G 

attains its minimum on G^ by compactness, and by the previous argument it is the global 
minimum on □ 
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Let iijsi denote the equidistributed probability on the set {ujj] 1 < j < N} of Lemma 3.6 
Lemma 3.7. The regions of the partition Vn of Q ean be indexed as in sueh way that 


Ch^ > M. \(jOj{0) — x\‘^dx. 


(18) 




Proof. Let C ^ collect the barycenters of the partition Vn, and let un denote the equidis¬ 
tributed probability on B]\f. One has 


TL2(Leb,eo #/iAr) < CeW2{ia^ i^n) = Cehjsf. 


VF 2 (^a/'? Leb) < CjyN 

Thus VL 2 (z>'Ar, eo #/iAr) < Cihjsi by Lemma [3^ This optimal transport problem between the 
discrete measures and eo#/iw determines an optimal assignment Ttv ^ B^^ represented 
by the indexation of Tw and Bjsf. Denoting by Pj G Vn the region of which hj is the 


barycenter we conclude that 

J \(jjj{0) — x\‘^dx = ^ J \bj — xf^dx+ W 2 {i^n^^o#Tn) ^ Ch% 


□ 


by 


^<j<N l<j<N^ 

For each 0 < i < T, let G N be the piecewise constant map on the partition Vn defined 

VI < j < Vx G Pj, mi{x) := ujj{i/T). 


Bound on the energy terms Hm^+i — m^||. Using Cauchy-Schwartz’s inequality we obtain 


T E ||m^+i - rui 

o<i<r 





\ujj{t)\^dt 


1 

N 



uj dpj{(jj) 


2 



f f \oj\‘^ dpj{(jj) dt = d^{s^^ s'"). 


Distance to incompressible maps. For any 1 < i < T, with t := i/T, one has 

inf \\mi - s|| = LF 2 (Leb,et = W 2 {et#p,et#pN) < CeW 2 {p,PN) = CehN> 


Boundary conditions. We make the assumption that s^ = Id, up to a minor modification of 
Lemma 3.7 (replace x with 5*(x) in ([T^). Lemma 3.7 then precisely states that \\mo — s*|p < 
Chj^, and the generalized boundary condition of Q states that uj{l) = s*((x;(0)) for /i- almost 
every cj G F. Denoting by Cq the Lipschitz regularity constant of 5* we obtain for any 1 < j < N 
and X G X 

|6c;j(l) — s*(x)p = J s''{(jj{0)) — s''{x) dpj{uj) < J |s*(6c;(0)) — s*(x)p (ipj(6L;) 

< Cq J — x\‘^ dpj{uj) < 2 Cq dpj{uj)-\-\ujj{0) — x\ 
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Therefore 


\mT — s 


* 112 


= [ \^j{'^)-s*{x)fdx 

s^c’o E 4 /;'-"" ' 

^<j<N \ 


— / \(jjj{0) — cj{0)f dpj{cj) + / \(jjj{0) — xf dx 

^<j<N \ / 

< 2Cl{ClW2{iJiN, Ij) + 11^0 - ■s*|p) < Ch%. 


Summation and final estimate. The value £{T^ N, A) of the minimum Q is 




||mj+i—mjll'^+A I ||mo — s*||"‘+ ||m,r — s*|P + inf ||mj — s|M < s*)+0(Th^A) 


0<i<T 


l<i<T 


4 Numerical experiments 

4.1 Minimization algorithm and choice of penalization 

We rely on a quasi-Newton method to compute a (local) minimum of the discretized problem 
(§. This means that we need to compute the value of the functional 

m G T ^ + AM|mo - + ||mT - s*|p + ^ d|(m^)y (19) 

o<z<r ^ i<i<T ^ 

and its gradient, where dg{m) = inf^^g ||m — sjp. The only difficulty is to evaluate the squared 
distance d| to the set of measure-preserving vector fields and its gradient. As explained in the 
introduction, Brenier’s Polar Factorization Theorem implies that for any vector valued function 
m G M, 

d|(m) = VF|(m#Leb,Leb). 

When m belongs to Mw, the measure m#Leb is finitely supported, and the computation of the 
Wasserstein distance can be performed using a semi-discrete optimal transport solver |AHA98[ 
IMerllf ILevl4| . The next proposition gives an explicit formulation for the gradient in term of 
the optimal transport plan. Recall that Mw is the set of piecewise constant functions on the 
tessellation Vn •— (^j)i<j<A^ of The diagonal Dw in Mw is the set of functions m in Mw 
such that m{Pj) = m{Pk) for some j ^ k. The set Mw \ is a dense open set in Mw- 

Proposition 4.1. The functional d^ is differentiable almost everywhere onMw and continuously 
differentiable on M^y \ The gradient of d^ at m ^ M^y \ ©w is explicit: with Xj = m{Pj), 

Vdl{m)\p_ = 2{xj - hary{T~^{xj))) ( 20 ) 

where T : X ^ m{X) is the piecewise constant optimal transport map between Leb and the 
finitely supported measure m#Leb and bary(S') = f^xdx/Leh{S) is the isobarycenter of S 

Proof. The functional T := — || • |p is concave as an infimum of linear functions: 

Tim) — dl(m) — IlmlP = inf \\m — slP — llmlp = inf \—2{m\s) + lIslPl , 

V / sv / II II II II II II L \ I / II II J ’ 
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where (•,•) denotes the L‘^{X) scalar product. This implies in particular that T and therefore 
dg is differentiable almost everywhere on Mw- Given m in Mw \ define xj = m{Pj) and 
let T : X ^ be the optimal transport plan from Leb to m# Leb = ^ ^xj • The 
transport plan is indeed always representable by a function when the source measure is absolutely 
continuous with respect to the Lebesgue measure. Let Vj — T~^{xj) be the partition of X induced 
by this transport plan. Then 


X{m) 


Leb, Leb) — \\m 



\Xj 


^dx 



\x\rdx 


where G{m) G Mw is the piecewise constant function on X given by G{m)\p. — —2baryCt^). 
For any m' in Mw and x'- — one has 


Leb, Leb) - 



= X{m) + {m! — m\G{m)) 


This shows that G{m) belongs to the superdifferential to T at m. In addition, by the continuity 
of optimal transport plans, the map m G Mw\®w G(m) is continuous. To summarize, on the 
open domain \®w the concave function X possesses a continuous selection of supergradient. 
This implies that X is of class on this domain, with VF{m) = G, and the result follows. □ 


Construction of the initial solution Since the discrete energy (19) is non-convex, the con¬ 
struction of the inital guess is important. We follow a time-refinement strategy already used by 
Brenier |Bre08| to construct a good initial guess. Assuming that we have already a local mini- 
mizer for T]^ — 2^ + use linear interpolation to construct an initial guess for Tfc+i = 2'=+i + l. 
The optimization is then performed from this inital guess, using a quasi-Newton algorithm for 
the energy (19). 


Choice of the penalization parameter The optimal choice of A in ( [T^ depends on the 
quantization dimension D — i?quant(/^) of fho generalized solution fi G Prob(f]) that one expects 
to recover: namely Aw = N~d ^ see the remark after ©• We call D the ffow dimension, and 
regard it as as the intrinsic dimensionality of the problem which determines its computational 
difficulty. For a classical solution, this dimension agrees with the ambient dimension i.e. D — d^ 
while for a non-deterministic solution the quantization dimension can be up to 2d. Intermediate 
dimensions d < D < 2d are also common |Bre89| . In our numerical experiments we set Aw = X 3, 
a decision justified a-posteriori by the numerical estimation of the quantization dimension of the 
computed solution, see Figure 

Note that the numerical error in © is governed (for a fixed number T of time step s) by 

the quantity A“^ -h and that /iw = 0{N~^) under the assumptions of Theorem 1.1 The 

- 

choice Xjsi — N oc thus yields a convergent scheme whenever a > d^ although convergence rates 
are improved if a is close to the ffow dimension T), so that Xn ^ Nd ^N- 
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4.2 Visualization of generalized solution 

The main interest of numerical experimentation is to visualize generalized solutions to Euler’s 
equation, or equivalently generalized geodesics between two measure-preserving diffeomorphisms 
5*, s* in S. 


4.2.1 Gradient of the pressure 

Consider a minimizer of the discretized energy Given i G {1,..., T — 1}, minimizes over 
Mw the functional m i-G T(||m — + ||^z+i — m|p) + \d?g{m). This gives 

- 2mi + rui+i) = T\Vd%{mi). (21) 


This equation is a discretized counterpart of the rule that the acceleration of a geodesic on an 
embedded manifold, is normal to that manifold (here S plays the role of the manifold, embedded 
in M, which is internally approximated by the linear space Mw)- The second order difference 
— 2mi + m^+i) approximates a second derivative in time. Comparing (21) to @ , we 
see that the right hand-side of (21) can be used as an estimation of (minus) the pressure of the 
gradient. 


4.2.2 Geometric data analysis 

As in the proof of Theorem [m the discrete minimizer of ( [T^ can converted to a collection of 
N piecewise-linear curves {cji,..., = r^v. We recall that the domain X is partitioned into 

N subdomains {Pj)i<j<N with equal area and we let ujj{ilT) G be the point corresponding 
to the restriction of rrii to the subdomain Pj, for each 0 < i < T. Figure illustrates this 
construction. We regard r^v as embedded in the Hilbert space IH := P^([0, 1],!^) which plays a 


natural role in the problem of interest, as in 
data analysis. 


and apply techniques from the field of geometric 


Clustering In order to better visualize the solution, we use the fc-means algorithm to divide 
the set rTV into k. A distinct particle color is attached to each cluster, see for instance Figure 
The /c-means algorithm consists in finding a local minimizer of the optimal quantization problem 


min — min 

^ i<i<k 
C(;Gr iv 


1*^ “ ^iWm 


( 22 ) 


using a simple fixed point algorithm, and to divide Tw into clusters {Ci)i<i<k with 


Ci = <cu e Tat; IIcc; — A||e = 


arg mm 
i<j<k 


\\(jj — t 


die 


, Ik automatically belong to Span(rw), hence to the d{T -\- l)-dimensional linear 
consisting of piecewise linear paths with nodes G at times t = i/T, 


Note that /i, • 
subspace of I 
0 < i <T. This makes (|22|) tractable. 


Box dimension A natural objective is to estimate the quantization dimension Pquant()n) of 
the generalized fiow ji G Prob(fl) minimizing the relaxed problem a. The probability measure 
/iw equidistributed on the set Tw approximates /i, see Proposition |1.2| hence we can expect 
the set Pw to also approximate supp(/i). The quantization dimension Pquant()n) is difficult to 
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estimate, but by Proposition 3.2 it admits the simpler upper bound L)box(supp(/i)). We estimate 


the latter by applying the furthest point sampling algorithm to the finite metric space Pat, which 
defines an ordering on the elements of Pat as follows: let 71 be an arbitrary point of Pat and 
define by induction 

7j+i := arg max d{j, {71,, 7*}) ( 23 ) 


As in Definition 3.1, denote by = r^(r n) is the smallest r > 0 such that Pat can be covered by 


i balls of radius r. For 1 <C i <C A^, the ratio log(i)/log(l/ri(rAr)) is expected to approximate 
log(i)/log(l/r^(supp(/i))) and thus the desired Dbox(supp/i). 


Lemma 4.2. Let £i := max^^p {ti, • • • ? 7 z}); where 7 ^ is defined as in (23). Then, 


( _ log(2) \ log(i) ^ log(0 ^ log(0 
V log(l/ej)/ log(l/£:j) “ log(l/ri) “ log(l/£:j) 


Proof. By construction, ri < £{. Moreover, the balls centered at the points 71 ,... , 7 ^ and with 
radius 6:^/2 are disjoint, so that > 6 :^/ 2 . □ 


4.3 Test cases and numerical results 

Our two testcases are constructed from two stationary solutions to Euler’s equation in 2D. Let 
s : M+ ^ S be a classical solution to Euler equation in Lagrangian coordinates starting from 
the identity map. We solve the discretized version ^ of the minimization problem ([^-([^, with 
5 * = 5 ( 0 ) = Id and 5 * := ^(tmax), where tmax > 0. For small values of tmax the solution to 
this boundary problem is simply the original classical flow s, but for larger values a completely 
different generalized flow is obtained. In this case the geodesic s in the space of the measure 
preserving diffeomorphisms is no longer the unique shortest path between its boundary values 
5* and 5*. The first classical behavior is guaranteed if the pressure hessian satisfies 

-< (vr/tmax)^ Id (24) 

uniformly on [0,tmax] x -A, see ^ and |Bre89| . In all the numerical experiments, the number of 
points is set to = 10 000 and the number of timesteps is T = 2^ + 1 = 17. 


4.3.1 Rotation of the disk 


On the unit disk D = {{xi^X 2 ) G x\ + x\ < 1}, the simplest stationary solution to Euler’s 
equation 0 is given by a time-independent pressure field and speed: 

p{xi,X2) = ]^{x\ + X 2 ), V{XI,X2) = {-X2, Xl). 


The corresponding Lagrangian flow s{t) is simply the rotation of angle t. The largest eigenvalue 
of is 1 at every point in D. Hence by (24) the flow of rotations is the unique minimizer to 
both the variational formulation ([^ and its relaxation ^ with boundary values s* = s(0) = Id 
and 5 * = ^(tmax), when t^ax < tt. Uniqueness is lost at the critical time t^ax = tt which 
corresponds to a rotation of angle tt, so that the final diffeomorphism becomes s* = s { 7 v ) = — Id. 
In this situation, the minimization problem ([^ has two classical solutions, namely the clockwise 
and counterclockwise rotations. The relaxation 0 has uncountably many generalized solutions 
such as, by linearity, superpositions of these two rotations. 

Another explicit example of generalized solution was discovered by Brenier |Bre89| : given a 
point X ^ D and a speed u, denote by ojx^v the curve ujx,v{t) — x cos(t) sin(t), t G [0,1]. Then, 
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Brenier’s solution is obtained as the pushforward by the map (x,u) ojx^v G ri of the measure 
on D xM? defined by 


li{dx^ dv) 


—T-L^{dx) 

TT 




n 


x\ 




where T-L^ denotes the /c-dimensional Hausdorff measure. In particular, the quantization dimen¬ 
sion of the solution is 3 = 2 -h 1. We refer to |BFS09| for more examples of optimal flows, and 
construct four dimensional one. Let be defined by combining (i) a classical rotation on the 
annulus D \ D{r), with D{r) = {x G \x\ < r} and (ii) Brenier’s solution rescaled by a factor 
r on the disc D{r). Then fij. is an optimal generalized flow of quantization dimension 3, whereas 
the averaged flow firdr is also optimal by linearity, and has quantization dimension 4. 


Numerical results The numerical solutions computed by our algorithm for the critical time 
^max = tt are highly non-deterministic. To see this, we select a small neighborhood around several 
points in the unit disk D and look at the trajectories emanating from this small neighborhood. 
As shown in Figure we can see that the trajectories emanating from each neighborhood fill 
up the disk. In addition, each indivual trajectory looks like an ellipse. Second, we estimate the 
box dimension of the support of the numerical solution (as explained in ^4.2). The estimated 
dimension is slightly above 3. 


4.3.2 Beltrami flow on the square 

On the unit square S = [—1/2,1/2]^, we consider the Beltrami flow constructed from the time- 
independent pressure and speed: 

p{xi,X2) = ^(sin(7rxi)^-h sin(7rx2)^) 

v{xi^X2) = (—cos(7rxi) sin(7rx2), sin(7rxi) cos(7rx2)) 

The maximum eigenvalue of is tt^, and |Bre89| implies that the associated flow is minimizing 
between = 5(0) = Id and 5* = 5(tmax) for tmax < 1- Because of the lack of symmetry, 
generalized solutions constructed from this flow are less understood than in the disk case. 


Numerical results Our numerical results suggest the following observations. First, as shown 
in Figure]^ the computed solutions with boundary values = Id and 5* = ^(tniax) approximate 
the classical flow if t^ax < 1, and are non-deterministic generalized flows if t^ax ^ L This 
suggests the sharpness of the bound given by |Bre89| . Interestingly, even for t > 1, the numerical 
solutions seem to remain deterministic in a neighborhood of the boundary of the cube. This can 
be seen more clearly in Figure where the particles have been divided into clusters using the 
/c-means algorithm (see ^4.2.2| ). 

The pressure gradient is estimated as in ^ 4.2.1 and is displayed in Figure These pic¬ 
tures seem to indicate a loss of regularity of the pressure near the initial and final times. This 
corroborates the result of |AF07] according to which the pressure beloners to Ll^i]0,n BV(X)). 

Figure [^suggests that the even for tmax = 1-5, the reconstructed solution for the Beltrami flow 
are more deterministic than the solution to the disk inversion. We estimate the box dimension 
of the support of the solution using the method explained in §4.2.2[ The result are displayed 


in Figure The estimated dimension is D = 2 for the deterministic solution (t^ax = 0.9) 
but it increases as the maximum time (and therefore the amount of non-determinism) increases. 
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Finally, we note that the estimated dimensions for ^ {1.1,1.3,1.5} seem to be strictly 
between 2 and 3, suggesting a fractal structure for the support of the solution. This would need 
to be confirmed by a mathematical study. 

Software. The software developed for generating the results presented in this article is publicly 
available at https://github.com/mrgt/EulerSemidiscrete 

Acknowledgement The authors thank Y. Brenier for constructive discussions and introducing 
them to the topic of Euler equations of inviscid incompressible fluids. 
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Figure 5: (First row) Beltrami flow in the unit square at various timesteps, a classical solution 
to Euler’s equation. The color of the particles depend on their initial position. (Second to fifth 
row) Generalized fluid flows that are reconstructed by our algorithm, using boundary conditions 
displayed in the first and last column. When < 1 we recover the classical flow, while for 
^max ^ 1 the solution is not classical any more and includes some mixing. 
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Figure 6: Using the fc-means cluster algorithm, we cluster the reconstructed trajectories for the 
Beltrami flow in the square with = 1-5 into 10 groups. This suggests that close to the 

boundary of the square the movement of particle is clockwise and deterministic while in the 
interior the movement is highly non-deterministic and counter-clockwise. 



(a) t = 0.0 


(b) t = 0.125t, 


(C) t = 0.25t„ax 




(d) t = 0.375t. 


(c) t — O.Stmax 


(f) t = 0.625t. 




Figure 7: Estimated pressure gradient for the Beltrami flow on the square with fmax = 1-5. 
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(a) {x,y) = (-0.7,0) 


(b) {x,y) = (-0.35,0) 


(c) {x,y) = (0,0) 



(d) {x, y) = (0.2,0) (e) {x, y) = (0.35, 0) (f) {x, y) = (0.5, 0) 


Figure 8: We select particles whose initial position lie in a small disk, and display their trajectories 
according to the computed solution to (Top) For the inversion of the unit disk (Bottom) 

For the Beltrami flow on the square, with tmax = 1-5. 




Figure 9: Estimation of the box counting dimension of the support of the computed solution, 
see i 4.2.2\ (Left) for the inversion of the unit disk (Right) Comparison between the estimated 
box counting dimensions of the solutions to the Beltrami flow on the square, depending on the 
maximum time. 
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